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A3 S TRACT 

A  queue  with  Poisson  arrivals  is  served  by 
exponential  servers.  When  a  server  becomes  free,  he  can 
serve  a  croup  of  customers  of  sice  ac  most  b  .  He  is  not 
allowed  to  process  a  group  of  sice  less  than  a  , 

The  actual  sice  of  each  croup  served  depends  on  the  number  of 
waiting  customers. 

The  queue  is  studied  as  a  Markov  process.  In  the  stable 
case,  the  stationary  probability  vector  of  this  process  has 
a  simple,  readily  computable  form.  Using  this  form,  the 
stationary  waiting  time  distribution,  which  is  analytically 
intractable,  may  be  expressed  in  an  algorithmically  useful 
form. 

Several  questions,  related  to  the  optimal  design  of  such 
a  service  system,  may  be  algorithmically  investigated. 

This  model  serves  to  illustrate  the  advantages  of  the  algorith¬ 
mic  approach.  The  design  criteria  are  not  analytically 
tractable  in  general. 

KEY  WORDS 

Multi-server  queues,  group  service,  optimal  design  of 
queues,  computational  probability. 
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Introduction 


The  purpose  of  this  paper  is  to  snow  -he  advantages  of 
ohe  algorithmic  approach  in  studying  the  optimal  design  problems 
for  a  simple,  but  useful  queueing  model.  Customers  arrive 
at  a  fixed  location  according  to  a  Poisson  process  of  rate 
•••  •  There  are  N  servers,  who  process  customers  in  groups. 

'•"hen  a  free  server  is  at  the  fixed  location  and  finds  i  _  a  >_  1 
customers  there,  he  removes  min  (i,b)  customers. 

The  parameter  b  denotes  the  maximum  allowable  size  of  a 
group  to  be  served  and  satisfies  b  >_  a  .  If  there  are  fewer 
than  a  customers  present,  all  free  servers  remain  at  the 
location  until  a  customers  are  present.  One  server  then 
leaves  with  a  group  of  a  customers. 

Services  are  thought  of  as  the  removal  of  groups  of 
customers  co  a  different  location.  The  times  between  the 
departures  and  subsequent  returns  of  a  server  are  assumed  to 
have  an  exponential  distribution  of  parameter  u  .  All 
service  times  are  mutually  independent  of  each  other  and  also 
of  the  interarrival  times  of  customers. 

There  are  a  number  of  practical  situations,  which  motivated 
the  study  of  this  queueing  model.  The  most  familiar  is  a 
taxi  stand  served  by  a  fleet  of  N  taxis.  The  customers 
arrive  at  the  stand  and  are  removed  in  groups  of  at  least  a 
and  at  most  b  .  The  quantity  1/p  is  then  the  mean  travel 


In  a  situation  of  greater  interest,  the-  N  servers  are 
locomotives  which  remove  wagon  loads  of  ore  from  a  mine  head 
to  ocher  locations.  The  paranecer  /.  is  now  the  rate  at 
which  wagon  loads  of  ore  (customers)  become  available  at 
the  station.  The  parameters  a  and  b  are  respectively  th 
minimum  and  maximum  allowable  numbers  of  wagons  per  train. 
One  may  also  view  the  N  servers  as  medical  evacuation 
helicopters  removing  casualties  from  a  field  hospiral  to  a 
more  remote  facility.  The  significance  of  each  of  the  five 
parameters  of  the  model  is  then  again  evident. 

The  design  problems  are  clear.  It  is  desirable  to 
choose  N  sufficiently  large  that  few  customers  will  be 
waiting  for  service,  yet  not  so  large  that  servers  spend 
an  inordinate  amount  of  time  waiting  idly  for  customers.  In 
applications ,  the  parameter  b  is  usually  fixed  by  techno¬ 
logical  constraints.  When  b  is  large,  it  may  be  advisable 
to  choose  the  parameter  a  considerably  smaller  than  b  ,  in 
order  to  expedite  service.  It  is  clearly  of  interest  to  be 
able  to  assess  the  influence  of  the  choice  of  a  on  the 
behavior  of  the  service  system.  Our  solution  procedure  is 
such  that  it  can  handle  problems  with  fairly  large  values  of 
N  and  b  ,  say  on  the  order  of  twenty  and  two  hundred 
respectively . 

This  queueing  model  may  be  denoted  by  the  symbols 
M/M (a ,b) /N  .  It  is  not  essential  that  service  times  be 


viewed  as  travel  times.  There  may  be  N  servers  who  remain 
fixed,  but  serve  customers  according  to  the  rule  specified 


by  the  threshold  values  a  and  b  .  Our  results  will  then 
provide  the  stationary  distribution  of  the  number  of  waiting 
customers  and  not  of  the  usual  queue  length.  There  are  a 
number  of  particular  cases  which  were  previously  treated. 

For  N  =  1  ,  the  model  is  discussed  for  general  service  time 
distributions,  dependent  on  the  batch  size,  ir.  i'euts  [1,2]. 

A  .general  procedure  to  analyze  the  queue  GI/?K(a,b)/l  is 
described  ir.  Chapter  4  of  [3j.  The  present  model  is  also  an 
instance  of  a  double-ended  queue.  It  adds  to  the  literature 
on  such  queues,  but  avoids  the  discouraginglv  formal  treat¬ 
ment,  which  hitherto  has  burdened  their  discussion. 

The  queue  is  studied  as  a  Markov  process  on  the  state 
space  E  =  [(i,j),  1  <  i  <  S  ,  C  j  <_  a-1;  U  >  -  •'  •  Th& 
process  is  in  the  state  (i,  j)  ,  1  £  i  <  N  ,  0  <_  j  a-1  , 
whenever  there  are  i  free  servers  and  only  j  customers. 

The  state  k  >_  0  ,  corresponds  to  the  case  where  all  servers 
are  occupied  and  k  customers  are  waiting.  It  is  convenient 
to  order  the  states  (i,j)  in  lexicographic  order  with  i  in 
descending  order.  The  states  k  _>  0  ,  are  listed  thereafter. 

The  structure  of  the  generator  Q  of  the  Markov  process 
is  important  to  the  analysis  of  the  model.  It  is  displayed 
for  N  =  4  ,  a  =  3  ,  and  b  =  6  .  The  general  structure  of 
Q  is  readily  apparent  from  this  particular  case. 

It  is  important  to  note  that  all  columns  labeled  k  >_  1  , 
contain  the  same  elements  i.e.  -  X  -  Ny  on  the  diagonal, 


immediately  above  and  Nv  in  the  row  labeled  k-^-b. 
The  upper  left  hand  submatrix  of  dimensions  (N+l)a  *  Ma 
has  the  form 
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where  the  matrices  C  and  E  are  .square  matrices  of  order 
a  ,  given  by 
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The  column,  labeled  0,  has  the  elements  -X-Ny  on  the 
diagonal,  X  immediately  above  and  the  entry  Nu  in  the 


rows  labeled  a,  a+1 


•  •  f 


b  . 


are  obtained 
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x..C  + 


lJ— N-l 


x  . 


E  +  X  .  :c- (N-j ) ul  J 


(N-j  +  l)ux.,  ,  =  0  , 
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for  2  <  j  <  N  . 


',4! 


x,  E  +  x1  [  C-  ( N-l )  ul  ]  +  xqNl  ( 1 , 
x,  a  ,a  -  "a+  N-  -  Nu$£(l-i ) ~1 (1 


5 

2,  X4e  +  Xn  (1-'  )  =  1  . 

j=l  J 


0  , 


This  system  has  a  (unique)  positive  solution. 


Proof 


The  modified  gecrr.etric  form  of  the  vector  x  and  the 
equilibrium  condition  (2)  follow  from  results  in  C.  B.  Winsten 
14],  also  discussed  in  greater  generality  in  M.  F.  Neuts  [3]. 

By  substitution  of  (3)  into  the  first  Na+1  steady-state 
equations  and  into  the  normalizing  equation,  we  obtain  the 
equations  (4).  It  is  readily  verified  that  the  coefficient 
matrix  of  the  first  Na+1  linear  equations  in  (4)  is  an 
irreducible  generator.  The  solution  vector  is  therefore 
determined  up  to  a  multiplicative  constant  and  all  its  components 
have  the  same  sign.  That  vector  is  now  uniquely  determined 
by  the  normalizing  equation. 


Corollary  1 


The  stationary  marcinal  density  ;o.’  of  the  number  o 

‘  j 

waiting  customers  is  given  by 

for  0  £  j  <_  a-1  , 

=  ,  for  '  >  a  . 

j  ~ 

The  stationary  probability  r^  that  j  servers  are 
free  is  given  by 

rj  =  x0  for  3  =  0 

=  x^e  ,  for  1  <  j  <  X. 


r 

-» 


The  case  N=1  may  be  solved  explicitly  in  terms  of 
Elementary  manipulations  lead  to 


for  0  <_  j  <_  a-1  , 


where 


x 


0 


For  a=l  ,  we  obtain  =  i-I  , 
which  corresponds  to  the  familiar  M/M/1  q 
on  computational  aspects  in  the  sequel  wi 


x-  =  Id  -  li  , 
ue.  All  comments 
be  reserved  for 


the  case  W"-2.  The  case  15=1  is  trivial. 


3.  The  Waitinc  Time  Distributions 


In  this  section,  we  consider  the  probability  W ■  x }  t 
j  <  b  ,  x  >_  0  ,  that  a  customer  arriving  to  the  stationary 

?s  service  before  or  at  time  x  and 
sup  of  size  j  .  Service  is  according 
served  discipline.  The  analytic 
functions  Wu  !  • )  are  exceedingly 

J 

ccm.c I icatad .  We  shall  shew  that,  in  contrast,  their  alccrithm.i 
characterization  is  straightforward. 

To  this  end,  we  first  consider  a  simple  absorption  time 
problem  ir.  a  .Markov  process  with  2b  states.  The  states  of 


queue  at  time  t=0 

,  en 

is  served  as  part 

of  a 

to  the  first-come, 

firs 

expressions  for  th 

e  it. as 

★  *  * 

this  process  are  denoted  by  l,2,...,b  and  1,2,  ...,  b 


*  ★ 

The  states  a  ,  . . . ,  b  are  absorbing.  The  rows  of  the 


*  * 

generator,  which  correspond  to  the  labels  a  ,  ...,  b  ,  are 


fore  identically  zero.  The  rows  with  labels  1,  2,  ...,  b 


1  ,  .  ..,  (a-1)  are  displayed  below,  again  for  the 

recresentative  case  a=3  ,  b=6  .  We  shall  consider  the 


orobabilitiss  i> .  .  *(x)  that  absorption  occurs  in  the  state 

i] 

★  ★  ★ 

j  €  {a  ,  ...,b  }  before  time  x  ,  given  that  the  process  is 


ec  in  zr.e  szaz* 


t  C  t  x.  /  ...f  ic;  ±  J 


Trie 


relation  of  these  probabilities  to  the  mass-functions  WM  { • ) 


will  soon  be  clear 


rl  N  — 


If  we  partition  the  (a-b-l)  >'  (2b)  array  into  a 

square  matrix  D  of  order  (a+fc-l)  and  a  matrix  D 

of  dimensions  (a-rb-1)  x  (b-a+1)  ,  then  it  is  readily  seen 

that  the  matrix  '  (x)  =  { f .  .* (x)  ;  is  civan  bv 

tj 

fX  -  - 
(5)  i (x)  =  I  exp (Du) D  du  =  Iexp(Dx)-IjD  D  , 

0 


for  x  >  0. 


Suppose  now  that  the  arriving  customer  finds  bv+r  ,  v  •_  1  , 

0  <_  r  <_  b-1  ,  customers  ahead  of  him  in  the  queue.  The  first 

bv  customers  will  be  removed  in  v  batches  of  size  b  . 

The  remaining  r  customers  and  the  arriving  customer  will 
be  served  together  in  a  batch  of  size  with  max  'a,  r-*-i)  <k<b . 

The  size  of  that  batch  will  depend  on  the  number  of  customers 
who  arrive  during  the  service  of  the  v  batches  and  of  those 

arriving  during  the  wait  for  a  free  server  .  The  following 
lemma  is  useful  in  the  sequel. 

Lemma  1 

The  probability  that  the  arriving  customer  finds  at  least 
b  customers  in  the  queue,  that  all  groups  of  size  b  ahead 
of  him  are  served  no  later  than  time  x  and  that  at  the 
beginning  of  the  service  of  the  last  such  group,  there  are 
j  ,  1  <  j  <  b  ,  customers  in  the  queue  who  will  be  served  in 
the  same  group  as  the  arriving  customer, 


is  given  by 


Ntdu 


j-1 


x 


Gj (x)  = 


i„x j-r-1 


_  H  -H  r*  n  —  *  n  f 

2,  x  ‘  i  expC-Ny  {l-5°)uj8  ^ - 

r=0  u  -1 0  (j-r-1)! 


(6) 


for  1  <  j  <  fc-1 


Gv  (x) 


,-x 

•  exo 
-  0 


Proof 


3y  the  law  of  total  probability,  we  see  that  for 
1  <  j  <  b-1  ,  the  probability,  described  in  the  statement 
of  the  lemma,  is  given  by 


sc 


(Nut) ' _1 
(v-l)l 


Nt 


This  readily  reduces  to  the  stated  expression.  The  proof 
for  j=b  ,  is  similar.  We  emphasize  that  the  arriving 
customer  is  now  included  in  the  count  described  by  the  index 

j  • 

We  can  now  describe  the  mass-functions  Wh  ( • )  for 
a  £  j  <  b  .  The  function  W  (•)  contains  an  extra  term, 

cL 

given  by 


Va(x) 


-  I 


N  a-2 

x.  .  +  I  I 
i=1  i , a-1  ii1  j£0 


x . 

l 


- A  u  (Xu 


a-j-2 


>  0 


(a-j-2)  ! 


X  du 


for  x  >_  0  .  This  term  corresponds  to  the  case  where  the 
arriving  customer  finds  one  or  more  free  servers. 


He  .nav  then  have  tc  wait  until  a  sufficient  t. urine r  of 
additional  customers  arrive  to  initiate  service  of  a  croup 
of  size  a  . 

The  states  {1,  b •  of  the  absorbing  Markov  process, 

which  was  defined  above,  correspond  to  the  situation  where 
all  complete  groups  of  size  b  ,  who  are  ahead  of  the  arriving 
customer,  have  been  served;  all  additional  arrivals  to  be 
included  in  the  next  group  have  been  counted  into  the  index 

j  ,  but  there  is  no  free  server  yet.  The  corresponding  state 

★ 

j  has  the  same  significance,  except  that  there  is  also  a 

free  server  present.  It  is  now  clear  that  absorption  into 

*  ★ 

one  of  the  states  ta  ,  ...,  b  }  signals  the  start  of  the 
service  of  a  group  of  the  corresponding  size.  This  group 
induces  the  arriving  customer. 

3y  e°  ,  we  denote  the  row  vector  (1,  C,  ...,  0)  of 
dimension  b-a+1  .  The  vector  0_  of  dimension  a+fc-1  is 
defined  by  ( 1 ,  i, ,  ...,  £  ^ ,  0,  ...,  0)  and  G(x) 

=  [G.(x),  ...,  G^fx),  0,  ...,  0],  The  row  vector  W(x)  of 
dimension  b-a+1  has  the  components  W  (x) ,  ...»  W.  (x)  . 

cl  D 

Theorem  2 

The  vector  W(x)  is  given  by 


(7)  W ( x )  =  V,  (x)  e  + 


X  _  ( X 

_9  exp  (Du)  D  du  +  I  G(u)  exp[D(x-u)] 
•'0  1 0 


for  x  >  0 


o  i 


I' 

I' 


Fro  or 


Formula  »/.'  is  imrediatf 


;.r  the  consideration  or 


several  absorption  tines.  The  first  term  corresponds  to 
-;.-5  case  are  the  arriving  customer  finds  one  or  mere  tree 


servers . 


adds  onlv  tc  W  ix,. 

*  a 


"he  second  tern  corresponds  to  the  rase 


arriving  customer  finds  j,  0  <_  "  b  ,  cos  toners  an.  ••c 
free  server  in  the  system. 

rp  j-'-jg  last  tern,  the  arriving  customer  fines  j  _  o 
customers  ahead  of  him.  The  waiting  tine  m  ...a.  ^ase  -s 
analysed  as  a  first  passage  time  in  the  preceding  dicussion 


Each  of  the  terms  in  Formula  (7)  may  be  computed  by 
the  numerical  solution  of  fairly  simple  systems  of  differential 
equations.  In  order  to  evaluate  V  <x)  .  we  rewrite  the 


defining  expression  as 
N  a-1 

V  lx)  =  I  J  Xij 
a  1=0  j =0  J 


«  1  -  x0(l-U 


t—7  . a-2- j 

-  \  U  (  '  U  )  i 

do  b  ;x*  TiWT 


exp  (Ax)  e 


for  x  >  0  . 


i 


The  ve::3 r  a  has  the  ccmcor.ents  a  -  7  >: .  ,  for 

J  i=l  13 

0  •_  2  i”-  •  The  matrix  of  order  a-1  is  the  same 

as  the  matrix  C  ,  de  -  i.nec  in  Section  1.  Clearly  V_  (x) 

a 

— '  » 

=  "I-",..  -  i-  >:•  e  ,  where  z^tx)  =  2q(x)-‘.  ,  for 


ht  '■•'her  to  evaluate  the  second  term,  we  write 


3\pOu  .  do  =  r '  exp  ■'  D:-:'  -  I 


It  is  now  clear  that  z^(x)  satisfies  the  system  of 
differential  ecuations 


t ,  ,  x)  =  z,  lx)  D  »  fa 


for  x  >  0  , 


=  u  .  .re  sevens  term 


•tia  \  /  is  tr.en 


ven  oy  z,  •  x;  0 


ixt  we  iniroou: 


ector-va i uec  function 


2 '■  - exp.D  x-u)  ’  du  , 


for  x  '■  0  . 


5y  multiplying  both  sides  by  exp(-Dx)  and  differentiating, 
we  see  that  the  vector  z_-,(x)  satisfies  the  system  of 
differential  equations 


i x )  =  z0  !x)D  G  (x) 


for  x  >  0  , 


with  ZjfO)  =  2  •  The  components  of  the  vector  G(x)  are 
evaluated  by  numerical  integration  in  the  expressions  (6). 


cf  the  vector  W(») 


,  a 


The  means  of  the  mass-f uncoions  W_.  ,  a  <  j  b  , 

nay  be  computed  from  Formula  {!)  .  Deno cine  larlace-Stielt',c“S 
transforms  by  the  addition  of  asterisks,  (7)  readily  yields 

*  -  *  _  7  - 

W  (s'  =  V  (s)e  +  -isZ-B’’  D  +  G  (s)  'sI-D)  ~D  , 


\  ,5.0 

a  — 


-  -  is:-D!'"D  -  G  ■.  s)  (Sl-D)  -G  <s'  ''Sl-D: 


!t  follows  that 


-  ?  ~  * 


(10)  W  =  -W  (0)  »  -V  ( 0 ) e  +  :r+G(»! :d  "D  G  !0,D  *D 


Clearlv 


1  ,  a-2 

-V*  (0)  =  f  2  'a-j-l'x. .  , 

’  i«l  j=0  i: 

I 

★ 

sc  that  it  suffices  to  evaluate  G  (0)  . 


r  rom  (6)  ,  we  obtain  routinely  that  for  1  <_  ; 


*  M’ :  J  r 

g,  (s>  =  ~  7  xn-; 


"i  —  7 

J r  _b+r 


-I  j-r 


s+X+Nu  (l-ib)J 


i  r  h  -p~r 

_  Nu  y  rbJ-i  7  i  X-t-Nu  (1-1  )  | 

^  r=0  j_s  +  X+Nu  ( 1-1^*)  i 


,  Nu  -b+j  j  i  U 


-1  1  r 


"X  x0 


r=l  s  +  A£ 


The  conditional  mean  waiting  time,  given  that  the  arriving 

customer  is  served  in  a  croup  of  size  j  is  given  bv  W-  W . ^  , 

3  3 

fcr  a  1  j  1  h  .  The  unconci  tic.nal  mean  waiting  time  is  given 


■  15)  W  e  =  -V*  (0)  -  ,  5  -r  G(®>  |  D  1e  -  G*  ( 0 )  e  . 

—  —  2  “  — 

Remark 


From  the  formulas  (12)  and  (13),  we  obtain  the  following 

alternate  extress ions  for  the  G.(x)  . 

3 

G ^  (x)  =  Xq lb  J  ^  Er  ( X ^  ^ ,  x )  ,  for  1  <_  j  <_  b-1  , 


Gb(x)»  ^f::0l2b  I  E  (XS"1.  *)*;U(x)+  ^  Ej_  (XC_1  (1~C)  ,  x) 
r=l  ! 


These  expressions  are  better  suited  for  numerical  computation 
tr.an  those  given  in  Formula  (£). 


From  Formula 
b 

(17)  G  .  (x)  =Nu 

j  =  l  3 

so  that 


(6) ,  we  obtain  by 


a  routine  summation  that 
expT -Nu  (1-^)  u]  du  , 

J  0 


(IS) 


If,  on  the  other  hand,  we  add  the  expressions  in  (14) ,  we 
obtain  after  some  calculations  that 


(19)  .  -G* 

The  expressions 
of  Equation  (1) 


(0 )  e 


in  (18)  and  (19)  are  equal,  since  by  virtue 
we  have 


>. 

This  provides  us  with 
calculations  involved 
mean  waiting  time. 


an  accuracy  check  or.  the  lengthy 
in  finding  the  expressions  for  the 


4.  Algorithmic  Aspects 

The  equilibrium  condition  (2)  readily  defines  the  smallest 
value  of  Nb  for  w7hich  the  queue  is  stable.  Equation  (1) 
is  rewritten  as 

b+T 

(20)  5  «  1  -  p  +  p€  , 

where  p  =  Nu(X+Nu)  ^  ,  and  is  solved  e.g.  by  Newrcr.'s  method. 
Defining  the  matrices  M(j)  ,  1  <_  j  <_  N  ,  by 

=  f (N-j)uI-C(f 1  , 


(21) 


M  ( j ) 


2 


(4)  may  be  rewritten  as 

x^J_1  e:  M(j)  ,  for  2<  j  £  N  , 

-r  x 2  e:  M(l)  . 

ive  performed  Gauss-Seidel  iteration  in  the  equations  (22), 
computed  the  next  value  of  Xq  from  the  penultimate  equation 
in  (4)  and  renormalized  each  solution  vector  by  the  last 
equation  in  (4).  Any  non-zero,  nonnegative  initial  vector 


Ix^,  ...,  x^,  Xg]  may  be  chosen  to  start 

the 

iterative 

solution. 

The  special  form  of  the 

matrices  C 

and 

E  may  b 

0 

used  to  economize  on  storage 

and  computer 

ional 

effort . 

On  ly 

the  first  component  of  each 

vector  Xj  Z 

dif  f 

ers  from 

zero 

and  it  requires  only  one  multiplication. 

The 

matrices 

(N-j)uI-C  ,  1  <  j  <  N  ,  are  of  the  form 


the  first  N  ecuations  in 


xx.=  M(N)  , 


(22)  = 
~3 


~  :  v»  _  4  _  1  >  . ,  v 

-  J  1  -*■  ' 


6 


-  A 


ere 


where  y,  =  ^  i  >./6)  J  ,  for  0  £  v  <_  a-1  .  It  is  hence  not 

necessary  -o  score  the  N  matrices  M(j)  .  The  computation 
of  the  stationary  probability  vector  of  the  matrix  Q  is 
both  stable  and  very  fast. 

The  vector  w (*■}  is  easily  computed  from  the  formulas 
(8)  and  (9) .  In  view  of  the  special  form  of  the  matrix  D  , 
the  vector  -[j?  +  Gl00)]  D  ^  is  obtained  with  very  little 
effort. 

As  was  discussed  in  Section  3,  the  vector  W(x)  of 
mass-f unctions  is  computed  by  the  numerical  integration  of 
simple  systems  of  linear  differential  equations  with  constant 
ccef f icients .  This  may  be  accomplished  by  any  one  of  a  number 
of  classical  methods.  The  special  structure  of  the  coefficient 
matrices  is  exploited  to  reduce  computation  time.  The  approach 
to  the  easily  computed  vector  W(°°)  provides  us  with  an 
accuracy  check  on  the  numerical  integration  procedure.  For 
examples  in  which  the  parameter  a  is  large  or  in  which  the 
queue  is  nearly  critical,  the  evaluation  of  W(x)  ,  x  >_  0  , 
requires  fairly  substantial  computation  times.  In  all  other 
cases,  the  approach  to  W(*)  is  very  rapid. 
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